/*
 * @Description: 
 * @Version: 2.0
 * @Autor: Seven
 * @Date: 2020-11-14 15:39:50
 * @LastEditors: Seven
 * @LastEditTime: 2020-11-28 15:06:13
 */
#include "FEMSpace.h"
#include <typeinfo>
#include <math.h>
#include <ctime>
#define pi 4.0*atan(1.0)
double f(double *p)
{
    return 2 * pi*  pi* sin(pi * p[0]) * sin(pi * p[1]);
}

double bc(double *p)
{
    return sin(pi * p[0]) * sin(pi * p[1]);
}

double f2(double *p)
{
    return 2 * pi*  pi* cos(pi * p[0]) * cos(pi * p[1]);
}

double bc2(double *p)
{
    return cos(pi * p[0]) * cos(pi * p[1]);
}
int main()
{
    clock_t starttime,endtime;
    starttime = clock();

    //g++ -o main test_possion.cpp -std=c++11 -I /usr/include/eigen3/
    RectangleDomain* r = new RectangleDomain({{0,0},{1,0},{1,1},{0,1}});
    int segmentx = 8;
    int segmenty = segmentx;
    //Mesh<2>* m = new Q1Mesh(r,{segmentx,segmenty});
    //Element<2>* e = new Quadrilateral_1_Element();
    Mesh<2>* m = new P1Mesh(r,{segmentx,segmenty});
    Element<2>* e = new Triangular_1_Element();
    //Mesh<2>* m = new Q2Mesh(r,{4,4});
    //Element<2>* e = new Quadrilateral_2_Element();
    //Mesh<2>* m = new P2Mesh(r,{8,8});
    //Element<2>* e = new Triangular_2_Element();
    Equation<2>* equ = new PossionEquation<2>();
    equ ->SetBoundaryConditionFunction(bc2);
    equ ->SetRightHandsTermFunction(f2);
    Possion_2D possionproblem(m,e,equ);
    possionproblem.AssembleStiffMatrix();
    possionproblem.AssembleRightHandsTerm();
    possionproblem.DealWithBoubdaryCondition();
    endtime = clock();
    std::cout << "方程离散用时："<<double(endtime - starttime) / CLOCKS_PER_SEC <<" s"  << std::endl;
    clock_t starttime2,endtime2;
    starttime2 = clock();
    possionproblem.Solve();
    endtime2 = clock();
    std::cout << "矩阵求解用时："<<double(endtime2 - starttime2) / CLOCKS_PER_SEC <<" s"  << std::endl;
    Eigen::VectorXd residual = possionproblem.A() * possionproblem.solution() - possionproblem.Rhs();
    double residual_l2 = residual.squaredNorm();
    double residual_max = residual.lpNorm<Eigen::Infinity>();
    std::cout<<"l2 residual is :" << residual_l2 << " ,residual max norm is :" << residual_max<<std::endl;
    Eigen::VectorXd Truesolution = Eigen::MatrixXd::Zero(m->n_point(),1);
    for(int j = 0; j<=segmenty;j++)
    {
        for(int i = 0;i <= segmentx;i++)
        {
            double point[2] = {i * m->x_h(),j * m->y_h()};
            Truesolution[i + (segmentx+1) * j] = bc2(point);
        }
    }
    Eigen::VectorXd error = Truesolution - possionproblem.solution();
    double error_l2 = error.norm() / segmentx;
    double error_max = error.lpNorm<Eigen::Infinity>();
    std::cout<<" Error l2 norm is :" << error_l2 << " ,error max norm is :" << error_max<<std::endl;
    //当网格为32×32时
    //方程离散用时：0.6122 s
    //矩阵求解用时：0.011739 s
    //l2 residual is :2.30088e-26 ,residual max norm is :2.24798e-14
    // l2 norm is :0.000401724 ,error max norm is :0.000803448
    
    //当网格为64×64时
    //方程离散用时：2.77136 s
    //矩阵求解用时：0.059209 s
    //l2 residual is :7.30099e-25 ,residual max norm is :9.02454e-14
    //Error l2 norm is :0.000100407 ,error max norm is :0.000200814
    
    //当网格为128×128时
    //方程离散用时：9.63577 s
    //矩阵求解用时：0.189469 s
    //l2 residual is :2.33393e-23 ,residual max norm is :3.61308e-13
    //Error l2 norm is :2.51002e-05 ,error max norm is :5.02004e-05

    //当网格为256×256时
    //方程离散用时：38.4457 s
    //矩阵求解用时：0.723362 s
    //l2 residual is :7.4755e-22 ,residual max norm is :1.44556e-12
    //Error l2 norm is :6.27496e-06 ,error max norm is :1.25499e-05

    //当网格为512×512时
    //方程离散用时：153.664 s
    //矩阵求解用时：2.88048 s
    //l2 residual is :2.39158e-20 ,residual max norm is :5.78256e-12
    //Error l2 norm is :1.56873e-06 ,error max norm is :3.13747e-06
    
    //当网格为1024×1024时
    //方程离散用时：747.213 s
    //矩阵求解用时：13.4154 s
    //l2 residual is :7.649e-19 ,residual max norm is :2.31306e-11
    //Error l2 norm is :3.92187e-07 ,error max norm is :7.84374e-07
}
